Spatial Data Classes in R

• Spatial entries in GIS can be represented as a vector or as a raster

• You see below how river can be represented as a vector (left) and as a raster (right)

Vector Layers

Vector layers are sets of geometries associated with non-spatial attributes

The geomteries are sequences of one or more coordinates, connected to form lines or polygons

The non-spatial attributes come as a table

Vector layers can be:

• points
• lines
• polygons

Spatial Data Structures in sp

Caveat

• We will avoid the use of sp, rgdal, and rgeos and focus on the sf package.

• The rgdal package is deprecated starting with 2023.

• You should avoid using those packages (i.e. you will get errors and maybe conflicts with other libraries).

• However, it is important to be aware of these packages, given how many forums, articles, and dependencies have been using them.

The sf package

sf is the standard package for working with vector data in R.

It relies on external software components: GDAL, GEOS, and PROJ.

It contains three classes:

• sfg — a single geometry
• sfc - a geometry column: a set of sfg geometries with CRS information
• sf - a layer which is an sfc column inside a data frame with non-spatial attributes.

The sf class

Here is a polygon layer with three features and six non-spatial attributes: BIR74, SID74, NWBIR74, BIR79, SID79, NWBIR79

The sf class

Mapping the layer

The sf class

Geometry

sfg is a single geometry which can be of multiple types

• st_point
• st_multipoint
• st_linestring
• st_multilinestring
• st_polygon
• st_multipolygon
• st_geometrycollection

Geometry

These are the seven most common Simple Feature geometry types

Geometry

These are the seven most common Simple Feature geometry types

Geometry

These are the seven most common Simple Feature geometry types

Geometry

These are the seven most common Simple Feature geometry types

Geometry

These are the seven most common Simple Feature geometry types

Geometry

These are the seven most common Simple Feature geometry types

Geometry

These are the seven most common Simple Feature geometry types

Geometry

The last type GEOMETRYCOLLECTION is rarely used

The Shapefile format does not support geometries of type GEOMETRYCOLLECTION

Geometry - Points

Let us create one point:

library(sf)
library(ggplot2)
pt1 = st_point(c(34.798443, 31.243288))

If we print the point, we get its WKT (Well-Known Text) representation

print(pt1)

POINT (34.79844 31.24329)

If we examine the class definition, we get:

class(pt1)

[1] "XY"    "POINT" "sfg"  

Geometry - Points

class(pt1)

[1] "XY"    "POINT" "sfg"  

• XY is the two-dimensional geometry of the point
• POINT is the geometry type (it could be POINT, MULTIPOINT, etc.)
• sfg is the general class (i.e. simple feature geometry)

Geometry - Points

This is what the point looks like:

plot(pt1)

Geometry - Polygon

This is how we create a polygon

a <- st_polygon(list(cbind(c(0,0,7.5,7.5,0),c(0,-1,-1,0,0))))

If we print the polygon, we get its WKT representation

print(a)

POLYGON ((0 0, 0 -1, 7.5 -1, 7.5 0, 0 0))

If we examine the class definition, we get:

class(a)

[1] "XY"      "POLYGON" "sfg"    

Geometry - Polygon

This is what that polygon looks like:

plot(a)

Geometry - Polygon 2

We can also make a more complicated polygon:

b = st_polygon(list(cbind(c(0,1,2,3,4,5,6,7,7,0),c(1,0,0.5,0,0,0.5,-0.5,-0.5,1,1))))
print(b)

POLYGON ((0 1, 1 0, 2 0.5, 3 0, 4 0, 5 0.5, 6 -0.5, 7 -0.5, 7 1, 0 1))

This is what the class looks like:

class(b)

[1] "XY"      "POLYGON" "sfg"    

Geometry - Polygon 2

This is what that polygon looks like:

plot(b)

Combining Geometries

We can combine geomteries using the c function

This will result into a single MULTIPOLYGON geometry, composed of all the shapes in the input

ab <- c(a, b)
print(ab)

MULTIPOLYGON (((0 0, 0 -1, 7.5 -1, 7.5 0, 0 0)), ((0 1, 1 0, 2 0.5, 3 0, 4 0, 5 0.5, 6 -0.5, 7 -0.5, 7 1, 0 1)))

If we check the class of ab, we see that it is a MULTIPOLYGON

class(ab)

[1] "XY"           "MULTIPOLYGON" "sfg"         

Combining Geometries

This is what the combined polygon looks like:

plot(ab)

Intersecting Geometries

A new geometry can be calculated by applying various functions to existing ones.

i <- st_intersection(a, b)
print(i)

GEOMETRYCOLLECTION (POLYGON ((7 0, 7 -0.5, 6 -0.5, 5.5 0, 7 0)), LINESTRING (4 0, 3 0), POINT (1 0))

Notice the difference from the combined geometries: ab

ab <- c(a, b)
print(ab)

MULTIPOLYGON (((0 0, 0 -1, 7.5 -1, 7.5 0, 0 0)), ((0 1, 1 0, 2 0.5, 3 0, 4 0, 5 0.5, 6 -0.5, 7 -0.5, 7 1, 0 1)))

Also, notice the difference between classes:

class(i)

[1] "XY"                 "GEOMETRYCOLLECTION" "sfg"               

class(ab)

[1] "XY"           "MULTIPOLYGON" "sfg"         

Intersecting Geometries

This is what the intersected geometries look like:

i <- st_intersection(a, b)
plot(i)

Intersecting Geometries

This is how they differ from the combined geometries:

ab <- c(a, b)
plot(ab)

Geometry column - sfc

Let us create two more points: pt2 and pt3.

pt2 = st_point(c(34.812831, 31.260284))
pt3 = st_point(c(35.011635, 31.068616))

Geometry objects can be collected into a geometry column - sfc, which is done with st_sfc

A geometry column also contains a Coordinate Reference System

Geometry column - sfc

We can combine pt1, pt2, and pt3.

Reminder: this is how we defined them:

pt1 = st_point(c(34.798443, 31.243288))
pt2 = st_point(c(34.812831, 31.260284))
pt3 = st_point(c(35.011635, 31.068616))

This is how we combine them:

combined <- st_sfc(pt1, pt2, pt3, crs = 4326)
combined

Geometry set for 3 features 
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 34.79844 ymin: 31.06862 xmax: 35.01164 ymax: 31.26028
Geodetic CRS:  WGS 84

POINT (34.79844 31.24329)

POINT (34.81283 31.26028)

POINT (35.01163 31.06862)

Note the we added 4326 - the code for the whole world.

Geometry column - sfc

The combined points look like below

plot(combined)

Layer sf

We can add attributes to the geometry in the following way:

name<-c("Point 1", "Point 2", "Point 3")
name_df<-data.frame(name)
name_df

     name
1 Point 1
2 Point 2
3 Point 3

Layer sf

We can combine the attribute table name_df and the geometry column geom (sfc) using st_sf, resulting in a layer

layer <- st_sf(name_df, combined)
layer

Simple feature collection with 3 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 34.79844 ymin: 31.06862 xmax: 35.01164 ymax: 31.26028
Geodetic CRS:  WGS 84
     name                  combined
1 Point 1 POINT (34.79844 31.24329)
2 Point 2 POINT (34.81283 31.26028)
3 Point 3 POINT (35.01163 31.06862)

Mapping the sf layer with plot

plot(layer)

Mapping the sf layer with ggplot

ggplot()+
  geom_sf(data=layer)+
  theme_bw()

sf properties: rows and columns

Some important properties of sf objects include rows and columns.

Number of rows:

nrow(df2)

[1] 2811

Number of columns:

ncol(df2)

[1] 3

We can obtain both rows and column by typing:

dim(df2)

[1] 2811    3

sf properties: bounding box coordinates

The st_bbox function returns the bounding box coordinates

st_bbox(df2)

    xmin     ymin     xmax     ymax 
12.21167 41.70574 12.77428 42.06974 

Reading vector layers

As discussed, the most common vector files are:

Reading vector layers

Here is how we would would read all these different files:

Reading a shape file

usa_cntry <- st_read(dsn="./data/usa_database/shape_files/gadm41_USA_0.shp", quiet = TRUE)

Reading vector layers

Reading a geojson file

usa_cntry = st_read('data/usa_database/geojson/gadm41_USA_0.geojson', quiet = TRUE)

Reading vector layers

Reading a from a geodatabase

usa_cntry <- st_read(dsn="./data/usa_database/usa_geodatabase.gdb",layer="gadm41_USA_0", quiet = TRUE)

Mapping the files

Conclusion

• Today we learned about vector layers: points, lines, polygons,
• We also learned about sources for country polygons: https://gadm.org
• We also learned about common vector files
• We learned about combining and intersecting geometries
• We also learned about EPSG codes
• We also learned about extracting layer components: geometry, attribute table, and coordinates.